import numpy as np

def gaussian_elimination(A, b):
    A = np.array(A, dtype=float)
    b = np.array(b, dtype=float)
    n = len(b)
    
    # 构造增广矩阵
    Ab = np.hstack([A, b.reshape(-1, 1)])
    
    for i in range(n):
        # 寻找主元行并交换
        max_row = np.argmax(np.abs(Ab[i:, i])) + i
        Ab[[i, max_row]] = Ab[[max_row, i]]
        
        # 对主元以下的行进行消元
        for j in range(i + 1, n):
            ratio = Ab[j, i] / Ab[i, i]
            Ab[j, i:] -= ratio * Ab[i, i:]
    
    # 回代求解
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - np.dot(Ab[i, i+1:n], x[i+1:n])) / Ab[i, i]
    
    return Ab, x

# 测试
A = [
    [31, -13, 0, 0, 0, -10, 0, 0, 0],
    [-13, 35, -9, 0, -11, 0, 0, 0, 0],
    [0, -9, 31, -10, 0, 0, 0, 0, 0],
    [0, 0, -10, 79, -30, 0, 0, 0, -9],
    [0, 0, 0, -30, 57, -7, 0, -5, 0],
    [0, 0, 0, 0, -7, 47, -30, 0, 0],
    [0, 0, 0, 0, 0, -30, 41, 0, 0],
    [0, 0, 0, 0, -5, 0, 0, 27, -2],
    [0, 0, 0, -9, 0, 0, 0, -2, 29]
]

b = [-15, 27, -23, 0, -20, 12, -7, 7, 10]

Ab, x = gaussian_elimination(A, b)
np.set_printoptions(precision=2, suppress=True)
print("消元后的增广矩阵：")
print(Ab)
print("方程的解 (列向量形式)：\n", x.reshape(-1, 1))

# 验证解的准确性
b_hat = np.dot(A, x)
assert np.allclose(b, b_hat, atol=1e-8), "解的验证失败：Ax ≠ b"

# 随机生成20阶方阵和非零向量
n = 20
A_random = np.random.randint(-100, 100, size=(n, n))
b_random = np.random.randint(-100, 100, size=n)
_, x_random = gaussian_elimination(A_random, b_random)
print("随机矩阵和向量求解的解 (列向量形式)：\n", x_random.reshape(-1, 1))

# 验证随机矩阵与向量的解
b_random_hat = np.dot(A_random, x_random)
assert np.allclose(b_random, b_random_hat, atol=1e-8), "随机解的验证失败：Ax ≠ b"